library(tidyverse)
library(broom) # tidy regression outputs
library(fixest) # fixest::feols() is very fast with big data and fixed effects
library(patchwork) # for combining plots
library(ggthemes) # for ggthemes
library(knitr) # for making regression tables in Latex
library(magrittr)
shift <- data.table::shift # for lagging/leading variables

mp_df <- readRDS("data/analysis/mp_df.rds")

#### Direct effect of protest on tweets ####

b1 <- feols(sum_ctweets_binary ~ fff_events_cumulative1 + sum_tweets |
              full_name + year_month,
            data = mp_df)
b2 <- feols(sum_ctweets_binary ~ fff_events_cumulative1 + sum_tweets + frontbench |
              full_name + year_month,
            data = mp_df)
b3 <- feols(sum_ctweets_binary ~ fff_events_cumulative1 + sum_tweets |
              full_name + year_week,
            data = mp_df)
b4 <- feols(sum_ctweets_binary ~ fff_events_cumulative1 + sum_tweets + frontbench |
              full_name + year_week,
            data = mp_df)

## Baseline effects table
baseline_table_binary <- bind_rows(b1 %>% tidy() %>% mutate(model = 1),
                                   b2 %>% tidy() %>% mutate(model = 2),
                                   b3 %>% tidy() %>% mutate(model = 3),
                                   b4 %>% tidy() %>% mutate(model = 4)) %>% 
  # Combine all model outputs
  set_colnames(c("term", "estimate", "std_error", "t_statistic", "p_value", "model")) %>% 
  # Retain only the protest variable
  filter(str_detect(term, "fff")) %>%
  # Denote specification/controls
  mutate(covariates = c("None", "Frontbench",
                        "None", "Frontbench")) %>% 
  # Extract model FEs
  mutate(fe1 = c(b1$fixef_vars[[1]],
                 b2$fixef_vars[[1]],
                 b3$fixef_vars[[1]],
                 b4$fixef_vars[[1]])) %>% 
  mutate(fe2 = c(b1$fixef_vars[[2]],
                 b2$fixef_vars[[2]],
                 b3$fixef_vars[[2]],
                 b4$fixef_vars[[2]])) %>% 
  # Extract model N
  mutate(nobs = c(b1$nobs,
                  b2$nobs,
                  b3$nobs,
                  b4$nobs)) %>%
  # Format numbers
  mutate_at(vars(estimate:std_error), ~format(round(., digits = 3), 
                                              nsmall = 3)) %>% 
  mutate_all(~as.character(.)) %>% 
  # Stack
  pivot_longer(names_to = "param", values_to = "value", 
               c(estimate:p_value, covariates:nobs)) %>% 
  # Paste parentheses around standard errors for table
  mutate(value = ifelse(param == "std_error", 
                        paste("(", value, ")", sep = ""), value)) %>% 
  # Clean up names and labels
  mutate(value = ifelse(param == "fe1" & value == "full_name", 
                        "MP", value)) %>%
  mutate(value = ifelse(param == "fe2" & value == "year_month", 
                        "Year--month", value)) %>% 
  mutate(value = ifelse(param == "fe2" & value == "year_week", 
                        "Year--week", value)) %>% 
  group_by(model) %>% 
  mutate(p_value = as.numeric(max(ifelse(param=="p_value", value, 0)))) %>% 
  mutate(value = ifelse(param == "estimate" & (p_value < 0.05 & p_value > 0.01),
                        paste(value, "*", sep = ""),
                        ifelse(param == "estimate" & (p_value < 0.01),
                               paste(value, "**", sep = ""),
                               value))) %>% 
  select(-p_value) %>% 
  ungroup() %>% 
  filter(param != "t_statistic",
         param != "p_value") %>% 
  # Spread to model columns
  pivot_wider(names_from = "model", values_from = "value", 
              names_prefix = "B") %>% 
  mutate(term = case_when(param == "estimate" ~ "FFF Protest",
                          param == "std_error" ~ "",
                          param == "covariates" ~ "Covariates",
                          param == "fe1" ~ "Unit fixed effect",
                          param == "fe2" ~ "time fixed effect",
                          param == "nobs" ~ "Observations")) %>% 
  select(-param) %>% 
  # Format as Latex table
  knitr::kable(., 
               format = "latex",
               booktabs = T,
               caption = "Main models: Binary measure of climate tweets",
               linesep = "")
# --> Copy baseline_table into Latex
baseline_table_binary

#### Comparing tweets and speeches, across different lag structures #### 

# Evaluate models with different lag structures for protests 

# First, in different windows: "cumulative" codings of the protest variable
# sw() argument: evaluates each of the variables in the brackets individually, so creates that many model in a list
fixest_1_binary <- feols(sum_ctweets_binary ~ sum_tweets + 
                    sw(fff_events_cumulative_lead2, fff_events_cumulative_lead1,
                       fff_events_cumulative0, 
                       fff_events_cumulative1, fff_events_cumulative2, 
                       fff_events_cumulative3, fff_events_cumulative4,
                       fff_events_cumulative5) |
                    full_name + year_week,
                  data = mp_df)

## Collect the model outputs from each step [sw()/csw()] in the fixest::feols models

sw_steps_binary <- seq_along(select(mp_df, starts_with("fff_events_cumulative")))

hold_fixest1_binary <- as.data.frame(matrix(NA))

## Bind all the fixest_1 model outputs together using loops
for (step in sw_steps_binary){
  int_fixest <- tidy(fixest_1_binary[[step]]) %>%   
    filter(str_detect(term, "fff")) %>% 
    mutate(step_number = step)
  hold_fixest1_binary <- bind_rows(hold_fixest1_binary, int_fixest)
}

# And tidy
hold_fixest1_binary <- hold_fixest1_binary %>% 
  select(-V1) %>% 
  slice(2:nrow(.)) %>% 
  set_colnames(c("term", "estimate", "std_error", "t_statistic", "p_value", "step_number")) %>% 
  mutate(days = as.integer(str_sub(term, start = -1L, end = -1L))) %>% 
  mutate(days = ifelse(str_detect(term, "lead"), -1*days, days)) %>% 
  mutate(model = "window")

## Speech models

## First, the most basic speech model to match the syntax of the Twitter models
fixest_i_binary <- feols(sum_cspchs_binary ~ sum_spchs +
                    sw(fff_events_cumulative_lead2, fff_events_cumulative_lead1,
                       fff_events_cumulative0, 
                       fff_events_cumulative1, fff_events_cumulative2, 
                       fff_events_cumulative3, fff_events_cumulative4,
                       fff_events_cumulative5) |
                    full_name + year_week,
                  data = mp_df)
# fixest_i

hold_fixest_i_binary <- as.data.frame(matrix(NA))

## Bind all the fixest_1 model outputs together
for (step in sw_steps_binary){
  int_fixest <- tidy(fixest_i_binary[[step]]) %>%   
    filter(str_detect(term, "fff")) %>% 
    mutate(step_number = step)
  hold_fixest_i_binary <- bind_rows(hold_fixest_i_binary, int_fixest)
}

# And tidy
hold_fixest_i_binary <- hold_fixest_i_binary %>% 
  select(-V1) %>% 
  slice(2:nrow(.)) %>% 
  set_colnames(c("term", "estimate", "std_error", "t_statistic", "p_value", "step_number")) %>% 
  mutate(days = as.integer(str_sub(term, start = -1L, end = -1L))) %>% 
  mutate(days = ifelse(str_detect(term, "lead"), -1*days, days)) %>% 
  mutate(model = "window")

# hold_fixest_i

## Plot with the Tweet models
figure_tweets_speeches_binary <- rbind(hold_fixest1_binary, hold_fixest_i_binary) %>% 
  mutate(ci_lower = estimate-1.96*std_error,
         ci_upper = estimate+1.96*std_error) %>% 
  mutate(outcome = rep(c("Tweets", "Speeches"), each = 8)) %>% 
  ggplot(., aes(days, estimate, color = outcome, group = outcome, shape = outcome)) +
  geom_hline(yintercept = 0, linetype = 2, color = "black") +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper),
                 position = position_dodge(width = 0.5), size=0.75) +
  geom_point(aes(shape = outcome), 
             position = position_dodge(width = 0.5), size=3) +
  scale_shape_manual(values = c(19,18)) +
  scale_x_continuous(breaks = seq(-2, 5, 1),
                     labels = c("-2\nW", "-1\nTh", 
                                "1\nFriday", "2\nF:Sa", 
                                "3\nF:Su", "4\nF:M", 
                                "5\nF:Tu",
                                "6\nF:W")) +
  scale_color_manual(values=c("#79b321","#344d0e")) +
  labs(x = "Days in local FFF protest window", 
       colour = "Outcome (binary)",
       shape = "Outcome (binary)",
       y = "Estimate and 95% confidence interval") +
  theme_tufte(base_family = "Arial") +
  theme(legend.position = "bottom",
        legend.text = element_text(size=15),
        legend.title = element_text(size=15),
        axis.title.x = element_text(size=15),
        axis.text.x = element_text(size=15),
        axis.title.y = element_text(size=15),
        axis.text.y = element_text(size=15),
        panel.border = element_rect(colour = "black", fill=NA, size=1),
        panel.grid.major = element_blank())
figure_tweets_speeches_binary
ggsave("plots/figa4.png", width = 8, height = 6, dpi = 600)

## Next, consider a set of models where time is treated differently
#' Want to be matching the opportunity structure that MPs have for speaking

# A first set of models as in the baseline model: tweets and speeches
fixest_speech_i_binary <- feols(sum_cspchs_binary ~ sum_spchs + frontbench +
                           fff_events_cumulative1 |
                           full_name + year_week,
                         data = mp_df)
fixest_speech_1_binary <- feols(sum_ctweets_binary ~ sum_tweets + frontbench +
                           fff_events_cumulative1 |
                           full_name + year_week,
                         data = mp_df)


## Second, consider a model where the daily speech data is collapsed to the week level

mp_df_week <- mp_df %>%
  select(date, year_month, year_week, parliament_sitting,
         full_name, lab1_con0, frontbench,
         sum_tweets, sum_ctweets_binary,
         sum_spchs, sum_cspchs_binary,
         fff_event) %>%
  # Create weeks starting on Friday to aggregate with
  mutate(week_friday = week(round_date(date, "week", week_start = 1)),
         year_week_friday = paste(year(date), week_friday, sep = "_")) %>%
  relocate(c(week_friday, year_week_friday), .after = year_week) %>% 
  group_by(full_name, year_week_friday) %>% 
  mutate(parliament_sitting_week = max(parliament_sitting, na.rm=T),
         week_sum_spchs = sum(sum_spchs, na.rm=T),
         week_sum_cspchs_binary = max(sum_cspchs_binary, na.rm=T), # sum() replaced with max()
         week_sum_tweets = sum(sum_tweets, na.rm=T),
         week_sum_ctweets_binary = max(sum_ctweets_binary, na.rm=T),
         week_fff_event = sum(fff_event),
         week_frontbench = max(frontbench, na.rm=T)) %>% 
  ungroup() %>% 
  # And aggregate to the weekly level
  select(year_month, year_week_friday, parliament_sitting_week, 
         full_name, lab1_con0, week_frontbench,
         starts_with("week")) %>% 
  ## NB. Should be about 85,000 observations, if not check that year_week and year_week_friday aren't both included with distinct()
  distinct() %>% 
  # And create a lagged protest variable
  group_by(full_name) %>% 
  mutate(week_fff_event_lag1 = shift(week_fff_event, type = "lag", n = 1)) %>% 
  # Remove rows with NAs in the treatment variable induced by lags/leads
  filter(!is.na(week_fff_event_lag1))

# Tweets at the week-level
fixest_speech_2_binary <- feols(week_sum_ctweets_binary ~ week_sum_tweets + week_frontbench + 
                           week_fff_event | # Don't need to lag the week because of how the weeks are defined
                           full_name + year_week_friday, 
                         data = mp_df_week)
# fixest_speech_2

# Speeches at the week-level
fixest_speech_ii_binary <- feols(week_sum_cspchs_binary ~ week_sum_spchs + week_frontbench +
                            week_fff_event |
                            full_name + year_week_friday,
                          data = mp_df_week)
# fixest_speech_ii

## Third, consider the day to the next parliamentary sitting instead of the week

mp_df_sitting <- mp_df %>% 
  # Group by the cumulative number of days that parliament has been sitting since session began
  group_by(tally_sitting) %>% 
  # Calculate characteristics of that sitting window
  mutate(min_window_year = min(year),
         min_window_month = min(month),
         min_window_week = min(week),
         window_length = as.numeric(max(ymd(date)) - min(ymd(date)))) %>% 
  ungroup() %>%  
  # Calculate whether an MP spoke/tweeted during that interval
  group_by(full_name, tally_sitting) %>% 
  mutate(sum_tweets_window = sum(sum_tweets),
         sum_ctweets_window_binary = max(sum_ctweets_binary), # sum() replaced with max()
         sum_spchs_window = sum(sum_spchs),
         sum_cspchs_window_binary = sum(sum_cspchs_binary)) %>%  # sum() replaced with max()
  # Calculate whether there were protests in that district during that interval
  mutate(sum_fff_events_window = sum(fff_event)) %>% 
  # Keep key variables
  select(date, parliament_sitting, tally_sitting,
         full_name, frontbench,
         starts_with("min_"), window_length, ends_with("_window"),
         sum_ctweets_window_binary, sum_cspchs_window_binary) %>% 
  ungroup() %>% 
  # Collapse data to the "parliament--sitting day period" --- a unit of opportunity
  distinct(tally_sitting, full_name, .keep_all=T) %>% 
  # Lag protest measure
  group_by(full_name) %>% 
  mutate(sum_fff_events_window_lag1 = shift(sum_fff_events_window, type = "lag", n = 1))

## Estimate the effect of climate protests on speaking in the next parliamentary sitting day

# Tweets at the sitting--day-level
fixest_speech_3_binary <- feols(sum_ctweets_window_binary ~ sum_tweets_window + 
                           frontbench + window_length +
                           sum_fff_events_window_lag1 |
                           full_name + min_window_year^min_window_week, # Combine year and week into a fixed effect
                         data = mp_df_sitting)
# fixest_speech_3

# Speech at the sitting--day-level
fixest_speech_iii_binary <- feols(sum_cspchs_window_binary ~ sum_spchs_window + 
                             frontbench + window_length +
                             sum_fff_events_window_lag1 |
                             full_name + min_window_year^min_window_week, # Combine year and week into a fixed effect
                           data = mp_df_sitting)
# fixest_speech_iii

## Put these together in a table
speech_table_binary <- bind_rows(fixest_speech_1_binary %>% tidy() %>% mutate(model = 5),
                          fixest_speech_i_binary %>% tidy() %>% mutate(model = 6),
                          fixest_speech_2_binary %>% tidy() %>% mutate(model = 7),
                          fixest_speech_ii_binary %>% tidy() %>% mutate(model = 8),
                          fixest_speech_3_binary %>% tidy() %>% mutate(model = 9), 
                          fixest_speech_iii_binary %>% tidy() %>% mutate(model = 10)) %>% 
  set_colnames(c("term", "estimate", "std_error", "t_statistic", "p_value", "model")) %>% 
  filter(str_detect(term, "fff")) %>%
  mutate(covariates = "Frontbench") %>% 
  mutate(outcome = rep(c("Tweets", "Speeches"), times = 3)) %>% 
  mutate(unit_observation = rep(c("Daily", "Weekly", "Sitting days"), 
                                each = 2)) %>% 
  mutate(fe1 = c(fixest_speech_1_binary$fixef_vars[[1]],
                 fixest_speech_i_binary$fixef_vars[[1]],
                 fixest_speech_2_binary$fixef_vars[[1]],
                 fixest_speech_ii_binary$fixef_vars[[1]],
                 fixest_speech_3_binary$fixef_vars[[1]],
                 fixest_speech_iii_binary$fixef_vars[[1]])) %>% 
  mutate(fe2 = c(fixest_speech_1_binary$fixef_vars[[2]],
                 fixest_speech_i_binary$fixef_vars[[2]],
                 fixest_speech_2_binary$fixef_vars[[2]],
                 fixest_speech_ii_binary$fixef_vars[[2]],
                 fixest_speech_3_binary$fixef_vars[[2]],
                 fixest_speech_iii_binary$fixef_vars[[2]])) %>% 
  mutate(nobs = c(fixest_speech_1_binary$nobs,
                  fixest_speech_i_binary$nobs,
                  fixest_speech_2_binary$nobs,
                  fixest_speech_ii_binary$nobs,
                  fixest_speech_3_binary$nobs,
                  fixest_speech_iii_binary$nobs)) %>% 
  mutate_at(vars(estimate:std_error), ~format(round(., digits = 3), 
                                              nsmall = 3)) %>% 
  mutate_all(~as.character(.)) %>% 
  pivot_longer(names_to = "param", values_to = "value", 
               c(estimate:p_value, covariates:nobs)) %>% 
  mutate(value = ifelse(param == "std_error", 
                        paste("(", value, ")", sep = ""), value)) %>% 
  mutate(value = ifelse(param == "fe1" & value == "full_name", 
                        "MP", value)) %>%
  mutate(value = ifelse(param == "fe2" & value == "year_week", 
                        "Year, week", value)) %>% 
  mutate(value = ifelse(param == "fe2" & value == "year_week_friday", 
                        "Year, week", value)) %>% 
  mutate(value = ifelse(param == "fe2" & value == "min_window_year^min_window_week", 
                        "Year, week", value)) %>% 
  group_by(model) %>% 
  mutate(p_value = as.numeric(max(ifelse(param=="p_value", value, 0)))) %>% 
  mutate(value = ifelse(param == "estimate" & (p_value < 0.05 & p_value > 0.01),
                        paste(value, "*", sep = ""),
                        ifelse(param == "estimate" & (p_value < 0.01),
                               paste(value, "**", sep = ""),
                               value))) %>% 
  select(-p_value) %>% 
  ungroup() %>% 
  filter(param != "t_statistic", 
         param != "p_value") %>% 
  mutate(term = "fff_protest") %>% 
  pivot_wider(names_from = "model", values_from = "value", 
              names_prefix = "B") %>% 
  mutate(term = case_when(param == "estimate" ~ "FFF Protest",
                          param == "std_error" ~ "",
                          param == "covariates" ~ "Covariates",
                          param == "outcome" ~ "Outcome",
                          param == "unit_observation" ~ "Unit of observation",
                          param == "fe1" ~ "Unit fixed effect",
                          param == "fe2" ~ "Time fixed effect",
                          param == "nobs" ~ "Observations")) %>% 
  select(-param) %>% 
  knitr::kable(., 
               format = "latex",
               booktabs = T,
               caption = "Online and offline speech: Binary measures of climate tweets and speeches",
               linesep = "")
# --> Copy speech_table into Latex
speech_table_binary


#### Heterogeneous effects ####

# Load in relevant MP-level covariates
mp_covariates <- read_csv("data/analysis/mp_hte_covariates.csv")

# Merge
mp_hte <- left_join(mp_df, mp_covariates, by = "full_name")


## Heterogeneous effects by party: Labour/Conservative

b_labcon_tweets <- feols(sum_ctweets_binary ~ sum_tweets + frontbench +
                             fff_events_cumulative1 * lab1_con0  |
                             full_name + year_week,
                           data = mp_df)

b_labcon_speeches <- feols(sum_cspchs_binary ~ sum_spchs + frontbench +
                               fff_events_cumulative1 * lab1_con0  |
                               full_name + year_week,
                             data = mp_df)

## Heterogeneous effects for the All-Parliament Group on Climate Change

b_apg_tweets <- feols(sum_ctweets_binary ~ sum_tweets + frontbench +
                          fff_events_cumulative1 * apg_climate_190731 | 
                          full_name + year_week,
                        data = mp_hte)

b_apg_speeches <- feols(sum_cspchs_binary ~ sum_spchs + frontbench +
                            fff_events_cumulative1 * apg_climate_190731 | 
                            full_name + year_week,
                          data = mp_hte)

## Heterogeneous effects for the Conservative Environment Network MPs v. Tories

b_cen_tweets <- feols(sum_ctweets_binary ~ sum_tweets + frontbench +
                          fff_events_cumulative1 * cen_caucus_191213  |
                          full_name + year_week,
                        data = mp_hte %>% filter(party_value == "Conservative"))

b_cen_speeches <- feols(sum_cspchs_binary ~ sum_spchs + frontbench +
                            fff_events_cumulative1 * cen_caucus_191213  |
                            full_name + year_week,
                          data = mp_hte %>% filter(party_value == "Conservative"))

## Heterogeneous effects for the Net Zero Scrutiny Group v. Tories

b_nzsg_tweets <- feols(sum_ctweets_binary ~ sum_tweets + frontbench +
                           fff_events_cumulative1 * netzero_scrutiny |
                           full_name + year_week,
                         data = mp_hte %>% filter(party_value == "Conservative"))

b_nzsg_speeches <- feols(sum_cspchs_binary ~ sum_spchs + frontbench +
                             fff_events_cumulative1 * netzero_scrutiny |
                             full_name + year_week,
                           data = mp_hte %>% filter(party_value == "Conservative"))

## Frontbench

b_frontbench_tweets <- feols(sum_ctweets_binary ~ sum_tweets + 
                                 frontbench * fff_events_cumulative1  |
                                 full_name + year_week,
                               data = filter(mp_df, !is.na(lab1_con0)))
b_frontbench_speeches <- feols(sum_cspchs_binary ~ sum_spchs + 
                                   frontbench * fff_events_cumulative1  |
                                   full_name + year_week,
                                 data = filter(mp_df, !is.na(lab1_con0)))


## Combine HTE tweets models into table
table_b_tweets <- bind_rows(b_apg_tweets %>% tidy() %>% mutate(model = 1),
                            b_labcon_tweets %>% tidy() %>% mutate(model = 2),
                            b_cen_tweets %>% tidy() %>% mutate(model = 3),
                            b_nzsg_tweets %>% tidy() %>% mutate(model = 4),
                            b_frontbench_tweets %>% tidy() %>% mutate(model = 5)) %>%
  set_colnames(c("term", "estimate", "std_error", "t_statistic", "p_value", "model")) %>%
  filter(str_detect(term, "fff")) %>%
  mutate_at(vars(estimate:std_error), ~format(round(., digits = 3),
                                              nsmall = 3)) %>%
  mutate_all(~as.character(.)) %>%
  pivot_longer(names_to = "param", values_to = "value",
               estimate:p_value) %>%
  mutate(value = ifelse(param == "std_error",
                        paste("(", value, ")", sep = ""), value)) %>%
  group_by(model, term) %>% 
  mutate(p_value = as.numeric(max(ifelse(param=="p_value", value, 0)))) %>% 
  mutate(value = ifelse(param == "estimate" & (p_value < 0.05 & p_value > 0.01),
                        paste(value, "*", sep = ""),
                        ifelse(param == "estimate" & (p_value < 0.01),
                               paste(value, "**", sep = ""),
                               value))) %>% 
  select(-p_value) %>% 
  ungroup() %>% 
  filter(param != "t_statistic",
         param != "p_value") %>%
  pivot_wider(names_from = "model", values_from = "value",
              names_prefix = "H") %>%
  add_row(term = "Fixed effect", param = "fe1",
          H1 = "MP",
          H2 = "MP",
          H3 = "MP",
          H4 = "MP",
          H5 = "MP") %>%
  add_row(term = "Fixed effect", param = "fe2",
          H1 = "Year, week",
          H2 = "Year, week",
          H3 = "Year, week",
          H4 = "Year, week",
          H5 = "Year, week") %>%
  add_row(term = "nobs", param = "nobs",
          H1 = as.character(b_apg_tweets$nobs),
          H2 = as.character(b_labcon_tweets$nobs),
          H3 = as.character(b_cen_tweets$nobs),
          H4 = as.character(b_nzsg_tweets$nobs),
          H5 = as.character(b_frontbench_tweets$nobs)) %>%
  mutate(term = case_when(term == "fff_events_cumulative1" ~ "FFF",
                          term == "fff_events_cumulative1:lab1_con0" ~ "FFF times Labour",
                          term == "fff_events_cumulative1:cen_caucus_191213" ~ "FFF times CEN",
                          term == "fff_events_cumulative1:netzero_scrutiny" ~ "FFF times NZSG",
                          term == "fff_events_cumulative1:apg_climate_190731" ~ "FFF times APG",
                          term == "frontbench:fff_events_cumulative1" ~ "FFF times Frontbench",
                          param == "fe1" ~ "Unit fixed effect",
                          param == "fe2" ~ "Time fixed effect",
                          param == "nobs" ~ "Observations")) %>%
  mutate(term = ifelse(param == "std_error", "", term)) %>% 
  select(-param) %>% 
  # select(-param, `APP-H1` = H1, `APP-H2` = H2, `APP-H3` = H3, `APP-H4` = H4) %>% 
  mutate_at(vars(2:6), ~ifelse(is.na(.), "", .)) %>% 
  {. ->> int_b_tweets } %>%
  # filter(term!="") %>% # Remove the standard error term
  knitr::kable(.,
               format = "simple",
               booktabs = T,
               linesep = "")

## Combine HTE speech models into table
table_b_speeches <- bind_rows(b_apg_speeches %>% tidy() %>% mutate(model = 6),
                              b_labcon_speeches %>% tidy() %>% mutate(model = 7),
                              b_cen_speeches %>% tidy() %>% mutate(model = 8),
                              b_nzsg_speeches %>% tidy() %>% mutate(model = 9),
                              b_frontbench_speeches %>% tidy() %>% mutate(model = 10)) %>%
  set_colnames(c("term", "estimate", "std_error", "t_statistic", "p_value", "model")) %>%
  filter(str_detect(term, "fff")) %>%
  mutate_at(vars(estimate:std_error), ~format(round(., digits = 3),
                                              nsmall = 3)) %>%
  mutate_all(~as.character(.)) %>%
  pivot_longer(names_to = "param", values_to = "value",
               estimate:p_value) %>%
  mutate(value = ifelse(param == "std_error",
                        paste("(", value, ")", sep = ""), value)) %>%
  group_by(model, term) %>% 
  mutate(p_value = as.numeric(max(ifelse(param=="p_value", value, 0)))) %>% 
  mutate(value = ifelse(param == "estimate" & (p_value < 0.05 & p_value > 0.01),
                        paste(value, "*", sep = ""),
                        ifelse(param == "estimate" & (p_value < 0.01),
                               paste(value, "**", sep = ""),
                               value))) %>% 
  select(-p_value) %>% 
  ungroup() %>% 
  filter(param != "t_statistic",
         param != "p_value") %>%
  pivot_wider(names_from = "model", values_from = "value",
              names_prefix = "H") %>%
  add_row(term = "Fixed effect", param = "fe1",
          H6 = "MP",
          H7 = "MP",
          H8 = "MP",
          H9 = "MP",
          H10 = "MP") %>%
  add_row(term = "Fixed effect", param = "fe2",
          H6 = "Year, week",
          H7 = "Year, week",
          H8 = "Year, week",
          H9 = "Year, week",
          H10 = "Year, week") %>%
  add_row(term = "nobs", param = "nobs",
          H6 = as.character(b_apg_speeches$nobs),
          H7 = as.character(b_labcon_speeches$nobs),
          H8 = as.character(b_cen_speeches$nobs),
          H9  = as.character(b_nzsg_speeches$nobs),
          H10  = as.character(b_frontbench_speeches$nobs)) %>%
  mutate(term = case_when(term == "fff_events_cumulative1" ~ "FFF",
                          term == "fff_events_cumulative1:lab1_con0" ~ "FFF times Labour",
                          term == "fff_events_cumulative1:cen_caucus_191213" ~ "FFF times CEN",
                          term == "fff_events_cumulative1:netzero_scrutiny" ~ "FFF times NZSG",
                          term == "fff_events_cumulative1:apg_climate_190731" ~ "FFF times APG",
                          term == "frontbench:fff_events_cumulative1" ~ "FFF times Frontbench",
                          param == "fe1" ~ "Unit fixed effect",
                          param == "fe2" ~ "Time fixed effect",
                          param == "nobs" ~ "Observations")) %>%
  mutate(term = ifelse(param == "std_error", "", term)) %>% 
  select(-param) %>% 
  mutate_at(vars(2:6), ~ifelse(is.na(.), "", .)) %>% 
  {. ->> int_b_speeches } %>%
  # filter(term!="") %>% # Remove the standard error term
  rbind(., c("Outcome", rep("Speeches", times = 4))) %>% 
  knitr::kable(.,
               format = "simple",
               booktabs = T,
               linesep = "")

## Combining tweets and speeches in one table
table_b_combined <- cbind(int_b_tweets[,1:2], int_b_speeches[,2],
                            int_b_tweets[,3], int_b_speeches[,3],
                            int_b_tweets[,4], int_b_speeches[,4],
                            int_b_tweets[,5], int_b_speeches[,5],
                            int_b_tweets[,6], int_b_speeches[,6]) %>% 
  rbind(., c("Outcome", rep(c("Tweets", "Speeches"), times = 4))) %>% 
  set_colnames(value = c("term", paste("B", seq(from = 13, to = 22, by = 1), sep = ""))) %>% 
  knitr::kable(.,
               format = "latex",
               booktabs = T,
               linesep = "")
table_b_combined

